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Abstract 

We compute the largest relaxation times for the totally asymmetric exclu- 
sion process (TASEP) with open boundary conditions with a DMRG method. 
This allows us to reach much larger system sizes than in previous numeri- 
cal studies. We are then able to show that the phenomenological theory 
of the domain wall indeed predicts correctly the largest relaxation time for 
large systems. Besides, we can obtain results even when the domain wall 
approach breaks down, and show that the KPZ dynamical exponent z = 3/2 
is recovered in the whole maximal current phase. 
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1 Introduction 

Several many particle systems are characterized by a steady mass transport. Ex- 
amples for this kind of systems can be found in biological transport [1] or vehicular 
traffic [2]. From a theoretical point of view these processes are of particular in- 
terest, because they exhibit generic non-equilibrium behavior. Due to the large 
number of important applications, many microscopic models for particle transport 
have been suggested in recent years [2-4]. 

Among these, the most important microscopic model for non-equilibrium par- 
ticle transport is the so called asymmetric exclusion process (ASEP) [5,6]. In this 
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model, particles jump on a one-dimensional lattice, either to the right (with prob- 
ability p dt) or to the left (with probability q dt), if the corresponding sites are 
empty. The model shows a number of generic effects [7] that are characteristic for 
non-equilibrium particle transport and maintain for the more specialized variants 
of the model [8, 9] . At the same time the ASEP is simple enough to obtain several 
exact results for the system, which is of great importance, because the general the- 
oretical framework of non-equilibrium physics is less developed. Exact results exist, 
e.g. for the stationary state of the system with periodic [6] and open boundary con- 
ditions [10-12]. The case of open boundary conditions is of special interest because 
one observes boundary induced phase transitions [7]. In this paper, we restrict 
ourselves to the totally asymmetric exclusion process (TASEP), that is q — 0. This 
case includes the most important phenomena but simplifies considerably the dis- 
cussion of the model. Open boundary conditions are implemented by two particle 
reservoirs that are coupled to the chain. The capacities of the reservoirs determine, 
together with the bulk hopping rates, the actual state of the system. 

The complexity of the open system is also reflected in the mathematical struc- 
ture of the stationary solution. While the steady state of the periodic system is 
given by a simple product measure [6] , it is highly non-trivial for the open system. 
Nevertheless, it can be calculated and it is possible to obtain several non-trivial 
quantities, e.g. current- or density fluctuations in the stationary state [10-12], or 
large deviation functions [13]. 

The dynamic properties of the open system are, however, more puzzling. Exact 
analytical results for the largest relaxation time r and the corresponding dynamical 
exponent z are so far only possible for the periodic chain by applying the Bethe 
ansatz [14, 15]. For the open chain estimates for r can be obtained from a phc- 
nomenological approach, that models directly the dynamics of the boundary layer 
separating the high and low density domains imposed by the particle reservoirs 
[16]. It has been shown, that this approach gives asymptotically correct results in 
a certain parameter regime [17]. For finite systems, as well as for general in- and 
output rates, relaxation times have to be calculated numerically. 

This has been done by U. Bilstein and B. Wehefritz [18], and later by M. 
Dudzihsky and G.M. Schiitz [19] who calculated the relaxation times by using exact 
diagonalization. This method is, however, restricted to very short chains (less than 
twenty), that are not in the asymptotic limit. In this work we used the density 
matrix renormalization group (DMRG) technique [20], that enables us to calculate 
the relaxation times for much larger system sizes, compared to [18, 19]. By treating 
large chains we have obtained more conclusive results for the dynamical exponent z 
in the maximal current phase and could also achieve convergence between the exact 
numerical values of the relaxation times and the estimates of the phenomenological 
approach. 

The paper is organized as follows. In the next section we discuss briefly the 
relevant physical concepts and the applied numerical techniques. In the third 
section we show the comparison of the domain wall predictions for finite systems 
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with our results. Section 4 is devoted to the special case of the disorder line 
a + P — 1. This section is followed by a discussion of the dynamical behavior when 
approaching the phase boundaries as well as in the maximum current phase. 

2 Analytic predictions and the DMRG method 
for non-equihbrium systems 

2.1 The TASEP with open boundary conditions 

For self-containedness we will repeat the definition of the model. The TASEP is 
defined on a one-dimensional lattice with L sites. The boundary sites of the chain 
are coupled to two particle reservoirs, one reservoir on the left that controls the 
particle input and a second on the right that governs the output of particles. 

We regard the process in continuous time (see [21] for a comparison of the 
different update procedures) , which corresponds to a random sequential update in 
computer simulations. If a link between sites i and i + 1 is selected, a particle 
located at i moves to site i + 1 if site ? + 1 is empty (for convenience we set the 
hopping rate to one). In case of choosing the link (0, 1) one introduces a particle 
with probability a if the first site is empty. Finally a particle may leave the system 
with probabihty /3, if the hnk (L, L + 1) is chosen. 

Each configuration a can be written in terms of boolean lattice gas variables 
(jj, i.e. (Tj = 0(1) if the site is empty (occupied). If we introduce an orthonormal 
basis |a) = |(Ti, (Tl) in the 2^-dimensional configuration space, we can define 
the probability vector \P{t)) as \P{t)) = "^^^j P{{o'},t)\a) . The time evolution of 
\P{t)) is determined by means of the master equation, that can be written as a 
Schrodinger equation in imaginary time [5, 21]: 

l\P{t))^-n\p{t)), (1) 

where 7i denotes the stochastic Hamiltonian. The matrix elements of Ti. arc the 
rates w{a — > a') for a transition a — > a'. Explicitly Tl is given by (cr|7i|cr') = 
—w{a' — > a) for the off-diagonal elements {a ^ a') and by 
{a\H\a) = J2{a}f^{a'} '^{'^ ~^ diagonal elements. 

The exact stationary solution of the master equation can, for general a, f3, be 
written as a product of infinite dimensional matrices [11]. This solution allowed to 
calculate the a, (3 dependence of the stationary quantities, e.g. the average flux. 
The results are summarized in the phase diagram of the system shown in fig. 1. 
Three phases can be distinguished by means of a different functional behavior of 
the fiow. In the low density phases AI and All, i.e. for a < f3, a < 1/2, the fiow 
J = J (a) is given by J = a{l — a) and analogously in the high density phases BI 
and BII by J = /?(1-/?). If both a > 1/2 and /5 > 1/2 (phase C) the capacity of the 
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Figure 1: Phase diagram of TASEP with random sequential update, depending 
on the input rate a and output rate /5. The insets show typical density profiles. 
Phase transitions are indicated by the solid lines. Along the diagonal dashed line 
a + /3 = 1, corresponding to the disorder line of the model, the density profile is 
flat. Finally, the transitions between the subphases AI / All (respectively BI / 
BII) are marked by dashed lines. 

particle reservoirs exceeds the capacity of the chain. Then the flux is independent 
of a and /3 and given by J = 1/4. This phase is called maximal current phase. 

Both the high and low density phases are divided into two subphases. In phase 
AI and BI, the capacities of both reservoirs are below the chain capacity. In phase 
All (BII), only the capacity of the exit (entrance) exceeds the chain capacity. 
This has, e.g. consequences for the asymptotics of the density profile [12], and, 
a question we address in this article, possibly also for the dynamics of the chain. 
Another important line is given for a + /? = 1. On this one-dimensional line the 
stationary solution is much simplified, i.e. it is given as product measure. 

Now we shall summarize known results for the relaxation times. For the periodic 
system, exact results for the energy-gap, i.e. the inverse of the largest relaxation 
time Ti, have been obtained by Bethe-ansatz techniques [14, 15]. The largest relax- 
ation time scales for any density of the system asymptotically as ri ~ L^, where 



As translational invariance is lost in the open system, the technique does not 
apply. Besides, the dynamics is profoundly modified by the presence of the two 
particle reservoirs, which impose the coexistence of two domains into the system. 



z = 3/2. 
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In case of the phases AI, BI, it is known that the two domains imposed by the 
reservoirs have a simple factorized structure. In this parameter regime the so- 
called domain wall (DW) theory can be applied [16]. The DW theory uses a coarse 
grained description of the dynamics of the process: Each particle reservoir that is 
coupled to the chain imposes independently a domain of a given constant density 
Pl and pr. The two domains are separated by a localized domain wall, which 
performs a biased random walk. The bias is due to the different capacity of the 
two reservoirs and can be calculated simply by using the conservation of mass. In 
a finite system the position of the domain wall is confined between two rcficcting 
walls. The description of the process allows to calculate the stationary as well as 
the fully time dependent probability distribution of the domain wall positions [19]. 
Then, it is in particular possible to estimate the largest relaxation times of the 
system, that are given by 



D+ + D- - 2VD+D-cos{7rn/{L + 1)) 



n 



1,2, 



(2) 



with 



/?(1-/?) 



a(l — a) 



l-a-l3' 1-a-p ^ ' 

These results are valid for a,^ < 0.5. The remaining parameter space has to 
be explored numerically. 



2.2 The DMRG method for stochastic models 

In the previous section we have shown that the master equation can be rewritten 
as a Schrodinger equation in imaginary time. Now, in order to calculate the longest 
relaxation time of the system, one has to calculate the two lowest eigenvalues of 
eq. (1) (the first eigenvalue is trivially zero, with an eigenvector corresponding 
to the stationary state). The eigenvalue calculations have to be done with very 
efficient diagonalization methods, in order to reach sufficiently large system sizes. 
This is possible by applying the DMRG method [22]. 

This method was first developed to study the properties of strongly correlated 
electrons [20]. It has recently been generalized in order to treat stochastic many 
particle systems [23,24]. 

The idea of the DMRG is the following. One starts with a small system (here 
12 sites) that can be diagonalized with standard numerical methods. Then one 
performs a large number of renormalization cycles in order to increase the system 
size. At each renormalization cycle, first the system is enlarged by adding 2 sites in 
the bulk of the system, i.e. in a place where it is less likely that the eigenmodes will 
be perturbed. The largest eigenvalues of the corresponding enlarged Hamiltonian 
are computed. These eigenvectors are used in order to construct the density matrix 
of the system, which will be used in the next stage. Second, in order to avoid an 
exponential growth of the hamiltonian, a projection onto the "most important" 



5 



modes has to be done. It has been shown that the best choice is to keep the 
m leading eigenvectors of the density matrix [20]. The parameter m has to be 
chosen small enough in order to allow a fast calculation of the eigenvectors, but 
large enough in order to obtain a high numerical precision. The error due to each 
truncation is well controlled [22]. 

Although no fundamental differences between the DMRG method for stochastic 
and quantum mechanical systems exist, one has to make a certain effort in order to 
overcome the numerical difficulties. This has been done following the suggestions 
of ref. [23,24]. We now discuss some details of our implementation. 

The main difference between the quantum mechanical and stochastic problem 
is that the stochastic hamiltonian is not hermitian. Therefore, the calculation of 
the eigenvectors of the hamiltonian required at each renormalization step cannot 
be done by standard diagonalization techniques (as the basic Lanczos or Davids- 
son algorithms) but one has to apply, e.g. the Arnoldi method or Lanczos for 
non-symmetric matrices. These methods are quite efficient but, compared to their 
analogues for symmetric matrices, less stable [25]. We use the Arnoldi method, 
that turned out to be numerically more stable than the Lanczos method for non- 
symmetric matrices. Nevertheless the method is numerically well controlled, be- 
cause the accuracy of the calculated eigenvector can be obtained via the residual 
norm, i.e. we check that it is indeed an eigenvector. Actually, we have used the 
Arnoldi method in an iterative way, where the initial guess is the result of the 
previous Arnoldi run, until the desired precision is reached. 

Apart from being non-hermitian, there is another difference between stochastic 
and quantum mechanical problems. It is well known that the lowest eigenvalue of 
stochastic hamiltonians is always zero, corresponding to the existence of a station- 
ary state, and that the corresponding left eigenvector |0;) has all its coordinates 
equal to 1. Therefore the first task will be to calculate the right eigenvector |0j.). As 
we are interested in the asymptotic dynamics of the system, we have to calculate 
the energy gap G(L) as well, i.e. the next smallest eigenvalue and the corresponding 
left and right eigenvectors and \lr). Then the longest relaxation time r = ri 
will be given by 

r = G{L)-\ (4) 

After the calculation of these eigenmodes, one has to eliminate the least impor- 
tant degrees of freedom of the system, in order to keep the system at a manageable 
size. It has been shown that the relevant degrees of freedom are the eigenmodes 
associated to the largest eigenvalues of the density matrix p. The latter is defined 
depending on how many eigenmodes we need to calculate. Here, for the calculation 
of the gap, we take 

P = lTr{\0i){0i\ + |0,)(0,| + + . (5) 

Note that it is always possible to construct a symmetric form for the density matrix. 
Therefore, standard algorithms can be applied in order to calculate all eigenvalues 
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and eigenvectors of p as for quantum systems. 

Another trick that allowed to gain accuracy on the gap calculation is to compute 
first the right eigenmode \0r) for the fundamental state (the left is known), and 
then to define a new hamiltonian 

H'(A)^H + A\0r){0i\. (6) 

This hamiltonian H' has, apart from the groundstate, the same spectrum as the 
original hamiltonian H. The gap is now given by the fundamental state of H' [24]. 

Note that one requirement of the DMRG method is that the requested eigen- 
vahic is well separated from the next one. This may not be true when the system 
size increases, and then the calculation becomes instable. The system size limita- 
tions come from these instabilities rather than from runtime or memory require- 
ments. So, the DMRG method either gives very accurate results for the eigenvalues, 
or does not converge at all. 

Finally we want to mention another particularity of the system. As rules are 
different at each end of the system (input or output), we cannot use the left/right 
symmetry of the system as people do for closed systems, and the left and right 
parts of the system have to be computed alternatively. 

3 Test of the Domain Wall Theory at Small Sys- 
tem Sizes 

The comparison between DW predictions and the exact analytical for the stationary 
state shows that they agree in the limit of large system sizes and a,P < 1/2 [17]. 
For finite systems, however, the coarse grained picture slightly deviates from the 
exact results. 

In case of the dynamical properties, so far no exact analytical results exist. 
Therefore the DW predictions have to be checked with numerical methods. 

In a previous work, we have compared time-dependent density profiles in a non- 
stationary regime, and found good agreement with the domain wall predictions 
[26]. Here, we calculate explicitly the relaxation times of the process, in order to 
compare with the predictions of the DW method. The same has been done already 
by using the non-symmetric Arnoldi method [18, 19], but the system sizes that one 
can treat by the non-symmetric Arnoldi method are too small (L < 16 in [19]) in 
order to obtain convergence with the DW predictions. Our aim is to improve this 
convergence thanks to the ability of DMRG to treat larger systems. 

First we consider cases with both a and /3 < 1/2, where the domain wall 
theory is expected to work well, though this was not fully verified by the exact 
numerical calculation, due to size limitations. Fig. 2 shows the comparison of 
the inverse relaxation time with the domain wall predictions in the high density 
phase {a — 0.2, /? = 0.1). For comparison we also included the result for L = 16 



7 



obtained in [19], which perfectly coincides with the DMRG calculations. Obviously 
the theoretical predictions and the numerical results are in good agreement, if the 
length of the chain reaches thirty sites. 
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Figure 2: Gap as a function of the system size, for a — 0.2 and (3 — 0.1. The 
solid line indicates the domain wall prediction, while symbols represent the DMRG 
results. The circle was obtained in [19] for L=16. Error bars for the DMRG results 
are not given, as they are much smaller than the size of the symbols. The inset 
shows in a log-log plot the difference between the two curves. 

For a = 0.4 and f3 = 0.2, the DW predictions for small system sizes are slightly 
less accurate, as shown in figure 3. This could be expected, as the density difference 
between the two domains Ap = 1 — /3 — a is smaller and therefore the width of the 
shock larger than in the previous case. If, however, the shock touches the bound- 
aries the approximation of microscopic states by ©-functions is less appropriate. 

For a or f3 larger than 0.5, the coexistence is not between a high- and low-density 
domain, but, e.g. for a > 1/2, between a maximal current and high-density domain. 
This case is much more complicated than the coexistence of domains in the phase 
AI or BI. It is known from the exact results, that long ranged correlations exist in 
the maximal current regime, which implies that the maximal current domain has 
no simple structure, and cannot be easily described. 

A naive attempt to describe the maximal current domain, which uses the fact 
that the flux is a constant J = 1/4 and would lead to exponentially decaying 
density profiles, is to build the density profile iteratively from the boundary with 
a mean field formula like j — pi{l — Pi+i). But this yields a wrong density profile. 
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Figure 3: Gap as a function of the system size, for a — 0.4 and /3 = 0.2. The 
sohd Une indicates the domain wall prediction, while symbols represent the DMRG 
results. The insets shows in a log-log plot the difference between the two curves. 

It was also suggested [19] to use the exact results of the density profile. However, 
this approach has intrinsic shortcomings. First the DW is no longer self-contained 
and, second, it is no longer possible to specify the relaxation times in a closed 
form, because the hopping rates of the random walk are now site-dependent. But 
above all, the description of the system states as the juxtaposition of two phases 
with an independent boundary layer between them breaks down. Indeed, the wall 
dynamics is coupled to the whole structure of the maximal current phase, due to 
the long ranged correlations. 

So we explore this parameter regime guided by two questions: (i) Do we obtain 
a finite relaxation time in the phases All jBll in the limit of large system sizes? 
(ii) And, if it is the case, is it possible to model the maximum current phase as a 
flat domain with an effective density Pe//? 

First we checked the domain wall predictions close to the transition line at 
a — 0.5, i.e. for a — 0.51 and 0.55 (see figure 4). In this parameter regime one 
already has long ranged correlations in the maximal current domain, but at the 
same time the magnitude of the deviations from a fiat domain is rather small. Our 
results indicate, that within this parameter regime the DW formulas (2-3) still 
lead to satisfying results for quite large system sizes. On the other hand, the DW 
prediction with a. ~ 1/2, which was conjectured at first to be the proper one for 
any a > 1/2 [19], clearly underestimates the gap for any system size. 
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For larger values of a (e.g. a — 0.65 in figure 4), the algebraic corrections of 
the density profile become non-negligible. Equations (2-3) do not give anymore a 
good estimate. 
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Figure 4: Gap as a function of the system size, for a — 0.51, 0.55, and 0.65, and 
13 — 0.3. The solid hues indicate the domain wall prediction, while the symbols 
represent the DMRG results. The domain wall prediction for a = 0.65 is the upper 
curve. The dotted line gives the domain wall prediction calculated for a = 1/2. 
The DMRG results ya{L) were interpolated using an algebraic fit function f{L) = 
aL^+c. The inset shows log[/(L) — c] as a function of log L (symbols) in comparison 
with \og[ya{L) — c] (sohd lines). 



Next, we checked whether our results are compatible with a finite relaxation 
time in the limit L ^ oo. Therefore we interpolated our results for G{L) using the 
functional form f{L) = aL^+c. We find a very good agreement between the fits and 
our numerical data, with a finite value of c (see the inset of fig. 4), which indicates 
that the relaxation times remain finite. Our estimates for the relaxation times are 
r(oo) = 116.(6) for a = 0.51, r(oo) = 101.(7) for a = 0.55, and r(oo) = 96.(1) for 
a = 0.65. 

A finite relaxation time for L — > oo implies at the same time, that the dynamical 
exponent vanishes. This can be checked systematically by calculating the size 
dependent dynamical exponent given as: 

In [G{L + 2)] - In [G{L)] 
- ln(L + 2) - ln(L) " 
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Figure 5: Size dependent dynamical exponents for different values oi a {(3 — 0.3). 
The results are compatible with z{oo) —0. 

Figure 5 shows that our results are compatible with z = for all values of a we 
took into account. So these results for a > 0.5 further support that the relaxation 
times are finite within these subphases. 

Now, one can simply take the values t{oo) in order to calculate an effective 
density of the maximal current domain. This procedure leads to DW predictions 
for small systems, which do not agree very well with our numerical findings. We 
thus believe that the non-trivial structure of the maximal current domain must be 
considered in modeling the domain wall motion. 

4 Disorder line 

For general a, (5 it has been shown that the stationary weights are given as products 
of infinite dimensional matrices. Nevertheless a one-dimensional line in the param- 
eter space exist, i.e. a + (3 = 1, where the stationary weights are simple products 
as for the periodic systems. These product states are in analogy with disordered 
states of quantum spin models. At the disordered hue the system is homogeneous, 
i.e. one can not identify two different domains. This implies that the relaxation at 
this particular line is not governed by the motion of the domain wall. 

Then it is an open question to know whether dynamical properties are changed 
on this line and how the dynamical properties compare to the periodic system. 
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For the periodic system, Bethe ansatz predicts the dynamical exponent z — 'i/2 
[14, 15]. Besides, the first non-zero eigenvalue has an imaginary part except for the 
density 1/2. The divergent relaxation times of the periodic system are, however, 
related to its translational invariance [27]. 

By contrast, for the open system it is expected that the relaxation times are 
finite if a 7^ 0.5, because the density fluctuations spread with a nonzero drift 
velocity [7]. 

First, we have applied the DMRG method to the special point a — (5 — 0.5, 
where the drift velocity of an excess density vanishes. 

1.50 I • ^ • ^ • ^ . ^ ■ 1 
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Figure 6: Dynamical exponent 2; as a function of the system size, ior a — (5 — 0.5. 

Fig. 6 shows that our numerical results for z{L), obtained up to L = 244, agree 
with the exponent z — 'i/2 predicted for the periodic system. 

We have also explored the remaining part of the disorder line. We present our 
simulation results in fig. 7. They seem to exclude the case z = 3/2, and indicate 
that the relaxation time would be finite in the large L limit {z{(yo) = 0). 

Although these results are consistent with the expected relaxation behavior, we 
cannot exclude that a branch corresponding to complex eigenvalues could cross our 
solution when L becomes large. The DMRG method does not allow to distinguish 
whether the numerical instabilities are due to the intersection or the convergence 
of two eigenvalue branches. 

Nevertheless, we state that our numerical results are in accordance with a finite 
relaxation time and, therefore, consistent with the physical picture which was pre- 
sented in [7] . Further support for this scenario comes from its position in the phase 
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Figure 7: Dynamical exponent z as a function of the system size, along the disorder 
line. The curves from bottom to top correspond to a = 0.1, 0.2, 0.3, 0.35, 0.4, and 
0.45. Error bars are of the same order as the line thickness. 

diagram. The disordered line touches a phase boundary only in a single point, i.e. 
ior a — f3 = 1/2. All other points of the disordered hne are located inside a phase 
where the relaxation times of the system are finite. Therefore, it would be rather 
counter-intuitive to observe a diverging relaxation time for these parameters. 

5 Phase transitions and the maximal current phase 

If the capacity of both particle reservoirs exceeds the capacity of the chain, the 
maximal current phase is realized. In the maximal current phase both localization 
lengths are infinite. This divergence of length scales leads, e.g. to an algebraic 
slope of the density profiles. 

In this phase, previous numerical calculations by Bilstein et al [18] for L < 20 
extrapolated to large systems seemed to indicate an exponent 3/2 in the whole 
maximal current phase. 

We have computed the dynamical exponent in various points in the maximal 
current phase. Results are presented in figure 8, for a constant (3 = 0.8, a being 
varied through the second order transition (a = 0.5) and inside the maximal current 
phase (we have checked that of course, results are the same when a and (3 are 
interchanged). We clearly see a transition when the maximal current phase is 
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Figure 8: Dynamical exponent 2; as a function of the system size, for /3 = 0.8. The 
left figure presents results obtained in the All phase, where the maximal current 
phase coexist with a high density phase. For a = 0.5, we are exactly on the 2nd 
order transition line. The right figure shows results within the maximal current 
phase. 



entered. On the transition line a — 0.5, we recover 2; = 3/2 in the large system size 
limit. Besides, from our numerical results, the infinite size dynamical exponent 
seems to be equal to 3/2 in the whole maximal current phase, confirming the 
extrapolation in [18]. 



6 Conclusion 

The use of a DMRG approach to compute the largest relaxation time in the TASEP 
model has allowed to give new strong evidence of the validity of the domain wall 
picture when reservoirs control the flow. 

For a or f3 larger than 1/2, the capacity of one of the reservoirs becomes larger 
than the chain capacity, and one of the coexisting domains is then a maximal 
current phase. Our numerical results indicate that the relaxation times t{oo) are 
finite for L —> 00 as in the phases AI/BI. As the bulk density in a maximal current 
phase is 1/2 independently of the boundaries, one could have expected that the 
relaxation time t{oo) would be independent of the most efficient reservoir (i.e. 
independent of (3 in All and a in BII). This is not the case, as boundary layers in 
the maximal current domain decay only algebraically, and lead to a site-dependence 
of the density. Besides, strong correlations that cannot be easily characterized exist 
in this phase. This means that the maximal current phase cannot be modeled as a 
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flat domain with density 1/2. 

But even if we model the maximal current domain as a domain of a constant 
effective density Peff, which leads to the estimated value of r(oo), we observe 
large deviations between DW predictions and relaxation times for finite system 
sizes - again due to the non trivial structure of the maximal current domain. The 
DW predictions would probably be much improved if one considers site dependent 
hopping rates, but they are difficult to obtain since the structure of an isolated 
maximal current domain is not known. 

Anyhow, as in phase AI and Bl, the domain wall theory still offers a simple 
explanation for a finite relaxation time of the system in phases All and BII. 

On the disorder line, this picture cannot be applied, as no domain wall can 
be identified anymore, and a qualitatively different dynamical behavior cannot be 
excluded. However, though a wall cannot be defined anymore, we could still trace 
the dynamics of the density fiuctuations by introducing a second class particle. 
The second class particle performs a biased motion to one of the boundaries for 
any density different from 1/2. Therefore one expects that the density fluctuations 
are also driven out of the system with a flnite velocity, which should lead to flnite 
relaxation times, if a 7^ 1/2 [7]. 

Contrary, in the periodic system, which has the same simple structure of sta- 
tionary state, the relaxation times diverge as L ~ L^^^ for arbitrary densities. The 
divergent relaxation times have been related to the translational invariance of the 
periodic system [27] and are, therefore, not expected in case of the open system, 
except if a ~ P — 1/2. 

Indeed, our numerical results support this scenario. Besides, we have shown 
that the dynamical exponent z = 3/2 is also recovered in the whole maximal 
current phase, for which no theoretical prediction exists. 

So we could summarize the expected behavior of the relaxation time as follows : 

(i) In the whole phases AI, All, BI, and BII, except on the transition line 
a = (3, our results are consistent with a converging relaxation time. Its limit value 
as L — s> 00 is known from the DW theory in AI and BI. 

(ii) On the line a — P, and in the whole maximal current phase, the relaxation 
time diverges when L — > 00. The associated dynamical exponent 2; is 2 on the line 
a = f3 and 3/2 in the maximal current phase and at the point a = f3 = 1/2. 

Further improvement could be obtained from the use of Finite Size algorithm for 
DMRG (FSM), as described by Carlon et al [23], in order to gain some precision 
on our calculations - and thus to reach larger system sizes - in the whole phase 
diagram, and especially near critical lines. 

It would also be interesting to apply the DMRG method to other models, which 
have a non-trivial but not necessarily site dependent domain structure, in order to 
check whether the DW theory still describes correctly the relaxation behavior of 
the system. 
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